load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;**************************************
begin

VAR_measure = getenv("VAR_measure")
ilevel = toint(getenv("ilevel"))

;VAR_measure = "VAR_stddev"
VAR_measure = "VAR_R95_R50"
ilevel = 4

;---------------------

Nmodel = 60
N1 = 35        ;CMIP5
N2 = 25        ;CMIP6

nlat = 180
nlon = 360

Nindex = 4
index_name = (/"R90","R95","R98","R99"/)

Nvar = 3
var_name = (/"pr1d","pr5d","pr30d"/)

Nlevel = 6
warming_levels = (/"0C","0.5degC","1degC","2degC","3degC","4degC"/)

Nseason = 4
seasons = (/"MAM","JJA","SON","DJF"/)


;--------------read change---------------
diri := "/WORK/zhangwx/EConstraint/data/"
change := new((/Nseason,Nvar,Nmodel,Nindex,nlat,nlon/),float)

do iseason = 0,Nseason-1
do ivar = 0,Nvar-1

  f1 := addfile(diri+"CMIP5/pr/"+var_name(ivar)+"/"+seasons(iseason)+"/RX_pct/multimodel/change_frequency_RX_pct_"+var_name(ivar)+"_multimodel_warming_levels_"+seasons(iseason)+".nc","r")
  change(iseason,ivar,0:N1-1,:,:,:) = f1->change_frq_RX_pct(:,ilevel,:,:,:)               ;(35model, 4index, lat, lon)

  f2 := addfile(diri+"CMIP6/pr/"+var_name(ivar)+"/"+seasons(iseason)+"/RX_pct/multimodel/change_frequency_RX_pct_"+var_name(ivar)+"_multimodel_warming_levels_"+seasons(iseason)+".nc","r")
  change(iseason,ivar,N1:N1+N2-1,:,:,:) = f2->change_frq_RX_pct(:,ilevel,:,:,:)        ;(25model 4index lat lon)

end do
end do

PR := change/100.+1             ;% change --> probability ratio
copy_VarCoords(change,PR)


;-------------read pd Variability-------------
variability := new((/Nseason,Nvar,Nmodel,nlat,nlon/),float)

do iseason = 0,Nseason-1
do ivar = 0,Nvar-1

  f1 := addfile(diri+"CMIP5/pr/"+var_name(ivar)+"/"+seasons(iseason)+"/variability/multimodel/variability_"+var_name(ivar)+"_multimodel_"+seasons(iseason)+"_1997-2014.nc","r")
  variability(iseason,ivar,0:N1-1,:,:) = f1->$VAR_measure$          ;(35model lat lon)

  f2 := addfile(diri+"CMIP6/pr/"+var_name(ivar)+"/"+seasons(iseason)+"/variability/multimodel/variability_"+var_name(ivar)+"_multimodel_"+seasons(iseason)+"_1997-2014.nc","r")
  variability(iseason,ivar,N1:N1+N2-1,:,:) = f2->$VAR_measure$       ;(25model lat lon)

end do
end do


;----------------read model warming periods-----------------
warming_periods := new(Nmodel,integer)

f1 := addfile(diri+"CMIP5/tas/GMST/multimodel/GMST_multimodel_warming_levels_year.nc","r")
warming_periods(0:N1-1) = f1->year_center(:,ilevel)                   ;(35model)

f2 := addfile(diri+"CMIP6/tas/GMST/multimodel/GMST_multimodel_warming_levels_year.nc","r")
warming_periods(N1:N1+N2-1) = f2->year_center(:,ilevel)               ;(25model)



;**************************************************
corr := new((/Nseason,Nvar,Nindex,nlat,nlon/),float)
corr = 0.

do iseason = 0,Nseason-1
do ivar = 0,Nvar-1
  do iindex = 0,Nindex-1
     corr(iseason,ivar,iindex,:,:) = spcorr_n( variability(iseason,ivar,:,:,:), PR(iseason,ivar,:,iindex,:,:), 0)
  end do
end do
end do


Nmodel_valid := num(.not.(ismissing(warming_periods)))          ;number of models reaching the warming level
prob := rtest(corr, Nmodel_valid, 0)

copy_VarCoords(change(:,:,0,:,:,:),corr)
copy_VarCoords(change(:,:,0,:,:,:),prob)

print(Nmodel_valid+"")



;***************************************************************
do ivar = 0,Nvar-1
do iindex = 1,1                      ;R95
;do iindex = 2,3                     ;R98 R99

wks := gsn_open_wks("pdf","/WORK/zhangwx/EConstraint/pic/correlation/onlyCMIP/frequency/Corr_map_Variability_PR_4season_"+warming_levels(ilevel)+"_"+VAR_measure+"_"+var_name(ivar)+"_"+index_name(iindex))

gsn_define_colormap(wks,"cmp_b2r")
gsn_reverse_colormap(wks)


plot := new((/Nseason/),graphic)
plot_mk := plot

res = True
res@gsnFrame = False
res@gsnDraw = False

res@cnFillOn = True
res@cnLinesOn = False
res@cnLineLabelsOn = False

;res@mpProjection = "Robinson"
;res@mpPerimOn = False

;res@gsnSpreadColors = True
;res@gsnSpreadColorStart = 63
;res@gsnSpreadColorEnd = 4
res@cnMissingValFillColor = "white"

res@tmXTOn = False
res@tmYROn = False

  ;res@mpMinLatF=-40
  ;res@mpMaxLatF=60
  res@mpCenterLonF=0
  ;res@mpShapeMode="FreeAspect"

  res@tmBorderThicknessF = 1.0
  res@tmXBMajorThicknessF = 1.0
  res@tmXBMinorThicknessF = 1.0
  res@tmYLMajorThicknessF = 1.0
  res@tmYLMinorThicknessF = 1.0

  res@tmXBLabelFontHeightF = .017
  res@tmYLLabelFontHeightF = .017

res@gsnLeftStringFontHeightF = 0.022
res@gsnRightStringFontHeightF = 0.022

;res@gsnStringFont = "helvetica-bold"
;res@gsnLeftStringParallelPosF = -0.06
;res@gsnLeftStringOrthogonalPosF = 0.15
;res@tiMainFont = "helvetica-bold"
;res@tiMainFontHeightF = 0.019
;res@tiMainOffsetYF = 0.003
;res@tiMainFontThicknessF = 0.5

;res@cnFillPalette = "BlueWhiteOrangeRed"
res@lbLabelBarOn = False

;res@tiMainFontHeightF = 0.023
;res@tiMainFontThicknessF = 2

res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels := (/-0.6,-0.4,-0.3,0,0.3,0.4,0.6/)


res@gsnLeftString = ""
res@gsnRightString = ""

leftst = (/"a","b","c","d","e","f"/)

do iseason = 0,Nseason-1
  res@gsnLeftString = seasons(iseason)
  plot(iseason) = gsn_csm_contour_map_ce(wks,corr(iseason,ivar,iindex,:,:),res)      ;(/Nseason,Nvar,Nindex,nlat,nlon/)
end do



;*********************
siglvl = 0.05

mkres = True

mkres@gsnDraw =False
mkres@gsnFrame = False
mkres@cnLinesOn = False
mkres@cnFillOn = True
mkres@cnLineLabelsOn = False
mkres@cnLevelSelectionMode = "ExplicitLevels"
mkres@cnMonoFillPattern = False
mkres@cnFillDotSizeF = 0.0012
mkres@cnInfoLabelOn = False
mkres@lbLabelBarOn = False
mkres@gsnLeftString = ""
mkres@gsnRightString = ""
mkres@cnLevels = (/siglvl/)
;mkres@cnFillColors = (/"transparent","maroon1"/)
mkres@cnFillColors = (/"black","transparent"/)
mkres@cnFillPatterns = (/17,-1/)
mkres@cnMonoFillScale = False
mkres@cnFillScales = (/0.8,1/)

do iseason = 0,Nseason-1
plot_mk(iseason) = gsn_csm_contour(wks,prob(iseason,ivar,iindex,:,:),mkres)
overlay(plot(iseason),plot_mk(iseason))
end do



;--------------------------
txres = True
txres@txFontHeightF = 0.017

;gsn_text_ndc(wks,var_name,fspan(0.21,0.85,3),fspan(0.84,0.84,3),txres)

txres@txAngleF = 90
;gsn_text_ndc(wks,index_name,fspan(0.03,0.03,4),fspan(0.74,0.26,4),txres)


;*******************************************
;regional
;*******************************************
lonW := (/0,0,80,115,0,260/)
lonE := (/360,360,150,180,25,300/)
latS := (/45,-70,45,30,50,50/)
latN := (/70,-45,70,45,70,70/)

subregions := (/"NH_extratropics","SH_extratropics","Northern_Asia","WNP_monsoon","Europe","Eastern_NAmerica"/)
subregions_plot := (/"NH mid-latitudes","SH mid-latitudes","Northern Asia","WNP monsoon region","Europe","Eastern North America"/)
domain := (/"ocean","ocean","land","ocean","land","land"/)
Nregion := dimsizes(lonW)

;--------------------------------------------------------
plres = True
plres@gsLineColor = "blue"
plres@gsLineDashPattern = 0
plres@gsLineThicknessF = 1.5

plot_pl := new((/Nregion,4/),graphic)
do iregion = 2,3                      ;JJA
  plot_pl(iregion,0) = gsn_add_polyline(wks,plot(1),(/lonW(iregion),lonE(iregion)/), (/latS(iregion),latS(iregion)/), plres )
  plot_pl(iregion,1) = gsn_add_polyline(wks,plot(1),(/lonW(iregion),lonE(iregion)/), (/latN(iregion),latN(iregion)/), plres )
  plot_pl(iregion,2) = gsn_add_polyline(wks,plot(1),(/lonW(iregion),lonW(iregion)/), (/latS(iregion),latN(iregion)/), plres )
  plot_pl(iregion,3) = gsn_add_polyline(wks,plot(1),(/lonE(iregion),lonE(iregion)/), (/latS(iregion),latN(iregion)/), plres )
end do

do iregion = 4,5                     ;DJF
  plot_pl(iregion,0) = gsn_add_polyline(wks,plot(3),(/lonW(iregion),lonE(iregion)/), (/latS(iregion),latS(iregion)/), plres )
  plot_pl(iregion,1) = gsn_add_polyline(wks,plot(3),(/lonW(iregion),lonE(iregion)/), (/latN(iregion),latN(iregion)/), plres )
  plot_pl(iregion,2) = gsn_add_polyline(wks,plot(3),(/lonW(iregion),lonW(iregion)/), (/latS(iregion),latN(iregion)/), plres )
  plot_pl(iregion,3) = gsn_add_polyline(wks,plot(3),(/lonE(iregion),lonE(iregion)/), (/latS(iregion),latN(iregion)/), plres )
end do



;------------
txres2 = True
txres2@txFont = "helvetica-bold"
txres2@txFontHeightF = 0.012
gsn_text_ndc(wks,leftst((/0,1/)),fspan(0.07,0.54,2),fspan(0.76,0.76,2),txres2)
gsn_text_ndc(wks,leftst((/2,3/)),fspan(0.07,0.54,2),fspan(0.50,0.50,2),txres2)


;--------------------
resP = True
resP@gsnPanelTop = 0.95
resP@gsnPanelLeft = 0.05
resP@gsnPanelLabelBar  = True
resP@gsnPanelYWhiteSpacePercent = 4.
resP@gsnPanelXWhiteSpacePercent = 4.
resP@lbLabelFontHeightF = 0.012
resP@pmLabelBarHeightF = 0.06


gsn_panel(wks,plot,(/2,2/),resP)


end do              ;{iindex}
end do              ;{ivar}


end
